First-principles theory of the orbital magnetization 
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Within density functional theory we compute the orbital magnetization for periodic systems 
evaluating a recently discovered Berry-phase formula. For the ferromagnetic metals Fe, Co, and Ni 
we explicitly calculate the contribution of the interstitial regions neglected so far in literature. We 
also use the orbital magnetization to compute the EPR g-tensor in paramagnetic systems. Here the 
new method can also be applied in cases where linear response (LR) theory fails, e.g. radicals and 
defects with an orbital-degenerate ground-state or those containing heavy atoms. 
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The electric polarization and the orbital magnetization 
are well known textbook topics in electromagnetism and 
solid state physics. While it is easy to compute their 
derivatives in an extended system, the electric polariza- 
tion and the orbital magnetization themselves are not 
easy to formulate in the thermodynamic limit, due to 
the unboundedness of the position operator. The prob- 
lem of the electric polarization has been solved in the 
'90s by the Modern Theory of Polarization (MTP) (HQ, 
which relates the electric polarization to the Berry phase 
of the electrons. A corresponding formula for the or- 
bital magnetization has been found very recently 0, Q 
showing that this genuine bulk quantity can be evaluated 
from the ground state Bloch wavefunctions of the peri- 
odic system. Since the discovery of the MTP, a wealth 
of papers have appeared reporting its successful appli- 
cations to first principles calculations of dielectric and 
piezoelectric properties 0|. On the other hand, ab-initio 
calculations of the orbital magnetization via the Berry 
phase formula have not been reported in literature yet, 
except than for simple tight-binding lattice models. 

The origin of the orbital magnetization in molecules 
and solids is time-reversal breaking caused by e.g. spin- 
orbit (SO) coupling. In ferromagnetic materials the or- 
bital magnetization is a not negligible contribution to the 
total magnetization. Several papers in literature 0, @] 
showed that the orbital magnetic moment of simple fer- 
romagnetic metals (Fe, Co and Ni) is strongly underesti- 
mated within density functional theory (DFT) if using 
the local density approximation (LDA) or generalized 
gradient approximations (GGA). Empirical corrections 
like the orbital polarization (OP) |7J have been thus em- 
ployed to obtain a better agreement with the experimen- 
tal values. Nevertheless it remains an interesting ques- 
tion if e.g. functionals beyond LDA/GGA would be able 
to describe the orbital magnetization correctly Q. All 
previous ab-initio calculations have been however carried 
out in the muffin tin (MT) approximation, i.e. comput- 
ing the orbital magnetization only in a spherical region 
centered on the atoms, neglecting the contribution of the 



interstitial region. 

In this letter, we present first principles DFT calcula- 
tions of the orbital magnetization by evaluating the re- 
cently discovered Berry phase formula 0, H|. For the 
ferromagnetic phases of Fe, Co, and Ni we show that the 
interstitial regions contribute by up to 50% to the orbital 
magnetic moments. So far neglected in the literature 
these contributions are thus shown to be one source for 
underestimated ab-initio values. Furthermore we make 
use of a relationship between the orbital magnetization 
and the electronic g-tensor that can be measured in elec- 
tron paramagnetic resonance (EPR) experiments Q . We 
propose a new non-perturbative method that is highly su- 
perior to existing linear response (LR) approaches [9l.ll0j|. 
since it can deal with systems in which spin-orbit cou- 
pling can not be described as a perturbation. 

The total (sum of spin and orbital) magnetization can 
be defined from the derivative energy E to t with respect 
to the magnetic field B 
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where /„ is the occupation of the eigenstate n and in the 
most general case the expectation value is to be taken on 
ground state spinors ip n . In the last equality we take ad- 
vantage of the Hellmann-Feynman theorem. The Hamil- 
tonian in atomic units is 



H = i[p + aA(r)] 2 + l/(r) + 



er • [W(r) x (p + aA(r))] 



(2) 



where we drop the trivial spin-Zeeman term, reducing the 
magnetization according Eq. ((T|) only to its orbital part. 
We use the symmetric gauge A(r) = |B x r. The last 
term in Eq. ([2J is the leading spin-orbit term, describing 
the on-site SO coupling (with fine structure constant a 
1 i and the abbreviation g' — 2(g e — 1) 0, [To|]) and cr 
are the Pauli matrices. We neglect the spin other orbit 



2 



(SOO) term, in general a small contribution to the orbital 
magnetization and to the g-tensor 
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By inserting Eq. @ in Eq. ((T|) we obtain: 



M 



/« (V"n I r x v I 4>n) , 



(3) 



where v = — z[r, with H and V 1 computed at B = 0. 
This expression can be directly evaluated in a finite sys- 
tem, but not in extended systems because of the un- 
boundedness of the position operator and of the contri- 
bution of itinerant surface currents |3j. However, in pe- 
riodic systems and in the thermodynamic limit, Eq. Q 
can rewritten as a bulk property 0, 0] : 
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k X 



(9kU n k\ X {H k + £nk - 2e F ) \dkUrik) (4) 



where Hk is the crystal Hamiltonian with B = 0, e„k and 
«„k are its eigenvalues and eigenvectors, eF is the Fermi 
level, N c is the number of cells in the system and N k the 
number of fc-points. 

Eq. ([3]) and Eq. Q are valid at an all-electron (AE) 
level. To compute the orbital magnetization within 
a pseudopotential (PS) approach, we recall that a PS 
Hamiltonian {%) reproduces by construction differences 
and derivatives of the total energy. Thus we can still ob- 
tain M, from Eq. ([1}, if we replace dTl/dTi and ip n by 
the corresponding PS quantities dH/dB and tp n . 

We obtain the PS Hamiltonian in presence of spin- 
orbit coupling and uniform magnetic field with the 
Gauge Including Projector Augmented Waves (GIPAW) 

where H is 



method [12J. In particular % = T^TLTb, 



given by Eq. ((2|) and Tb is the GIPAW transformation 
[Eq. (16) of Ref. [H. If the AE and PS partial waves 
have the same norm the GIPAW hamitonian is given by 



ft — ft t rtgo t ft t ft-so 



0(B 2 



where 



W (0) = ip 2 + ^ ps (r) + ^L 
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^SO = T a 



(W ps (r) x p) + 
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(5) 



(6) 



H (1) = fB^L + 5:Rxi[r,^ L ]) (7) 



H<£<> = f-a 3 B- rx(<rxvg+^ 



R 



(8) 



Here V ps and V^ L are the local part, and the non-local 
part in separable form of the norm-conserving PS. 
and E^ L are the separable non-local GIPAW projectors, 
accounting respectively for the so-called paramagnetic 
and diamagnetic contributions (l3j of the atomic site R. 



Inserting H W + ftgo in E q- © 

we obtain: 

M = Mbarc + AMbaro + AM para + AM dia (9) 
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AM dia = 



g'a 3 
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^^r)xi[r-R,C]) (11 
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where (...) stands for / nk ( 

^nk I • • • I ^nk) ■ 

In a periodic system Mbarc can be nicely calculated by 
evaluating Eq. (gj) for the GIPAW Hamiltonian H and 
corresponding PS eigenvectors u„k and eigenvalues e n k- 
All the reconstruction terms, Eqs. (ITlTtT3|) . can be easily 
evaluated in extended systems, since the non-local oper- 
ators Vj^ L , and E^ L act only inside finite spherical 
regions, centered around each atom. 

The approach presented so far allows the calculation of 
the orbital magnetization in a general PS scheme includ- 
ing non-collinear spin-polarization. In this work for the 
sake of simplicity we use a collinear implementation. All 
expectation values are evaluated by assuming decoupled 
spin channels along the spin direction e. In particular 
all the spinors are eigenvectors of er ■ e and the local 
and total spin (S = Se) are aligned along e. Since the 
choice of e changes the spin-orbit coupling, the orbital 
magnetization is a function of e. In ferromagnets, each 
spin-direction e is characterized by a corresponding total 
energy, whereby the minimum of the total energy with 
respect to e defines the preferred direction of the spin- 
alignment, the so-called easy axis of the ferromagnet. 

We implemented our method in the Quantum- 
Espresso plane wave code [14j. We use standard norm- 
conserving pseudopotentials [l5| with two GIPAW pro- 
jectors per angular momentum channel. Using spin po- 
larized LDA [16[ and PBE 17 1 functionals, we perform 
standard SCF calculations including the SO term of 
Eq. ((6]) in the collinear approximation within the Hamil- 
tonian. Then we evaluate the orbital magnetization, ac- 
cording to the Eqs. H9I11H13P and Eq. (g]) for M ba rc- 
We neglect any explicit dependence of the exchange- 
correlation functional on the current density. In practice, 
spin-current density-funntional theory (SCDFT) calcu- 
lations have shown to produce negligible corrections to 
the orbital magnetization Q. We compute M(e) with 
e along easy axis and along other selected directions. 



3 



TABLE I: Orbital magnetic moments M(e) -e in \xb per atom 
of ferromagnetic metals parallel to the spin, for different spin 
orientations e. * denotes the experimental easy axis. The 
interstitial contribution is defined by the difference between 
M(e) and M^ b T = E„ k L, <k(r)r x (-iV + k)u„ k (r)dr 
where £l a is a atom-centered sphere of radius _Rmt=2.0 rbohr- 
All theoretical values are based on the gradient corrected PBE 
functional. The decomposition of M according to Eq. @ is 
shown Tab. II of the auxiliary material. 



Metal 


e 


Expt. 


This method 


FLAPW 






[20J 


Total 


Interstitial 


MT 


[5] 


6cc-Fe 


[001]* 


0.081 


0.0658 


0.0225 


0.0433 


0.045 


&cc-Fe 


[111] 




0.0660 


0.0216 


0.0444 




fcc-Co 


[111]* 


0.120 


0.0756 


0.0122 


0.0634 


0.073 


fcc-Co 


[001] 




0.0660 


0.0064 


0.0596 




hep-Co 


[001]* 


0.133 


0.0957 


0.0089 


0.0868 




hep-Co 


[100] 




0.0867 


0.0068 


0.0799 




fcc-Ni 


[111]* 


0.053 


0.0519 


0.0008 


0.0511 


0.050 


/cc-Ni 


[001] 




0.0556 


0.0047 


0.0509 





The k-derivative of the Bloch wave functions can be ac- 
curately evaluated by either a covariant finite difference 
formula [l8[ or by the k ■ p method [3| . For insulating 
systems both methods provide exactly the same results; 
for metallic systems the covariant derivative is more in- 
volved and we apply just the k ■ p method. 

For the ferromagnetic Fe, Co and Ni, calculations are 
carried out at the experimental lattice constants. We 
consider 4s and Ad states in the valence with non-linear 
core-correction, use a relatively low cutoff of 90 Ry. In 
the case of Fe, the results do not change by more than 
1% by including 3 s and 3p in valence and working at 120 
Ry. We use a Marzari-Vanderbilt cold smearing of 0.01 
Ry. We carefully test our calculations for fc-point con- 
vergence. In all cases a 28x28x28 mesh yields converged 
results within ±0.0001/is. 

Tab. H] reports our results for the orbital magnetiza- 
tion of the three metals Fe, Co and Ni, together with 
experimental values and a recent calculation performed 
by FLAPW @. All theoretical data was obtained us- 
ing the PBE functional. LDA gives within ±0.003m b th e 
same values (see Tab. II in the additional material [21|). 
In order to evaluate the contribution of the interstitial 
regions neglected so far in the literature (as in [5, 6]), 
we have equally computed (a/2) (L) only inside atomic 
spheres. Except for fcc-Co our results agree very well 
with FLAPW calculations. For Ni the influence of con- 
tributions is indeed negligible, explaining the agreement 
of early DFT calculations in this case. For the other fer- 
romagnets however it becomes evident from Tab. U that 
these contributions can by no means be neglected. For 
Fe e.g. the interstitial contribution is about 50% of that 
inside a MT sphere and thus leads to considerably im- 
proved ab-initio values. This result indicates the impor- 
tance of the contributions from the interstitial regions 



TABLE II: Principal values Ag in ppm for the diatomic 
molecules of the RnF-family calculated by linear response 
(LR) [13] and with the current method. || is symmetry axis 
of the dimer. Ag(AM) gives the contributions of AMbarc, 
AM para , AMdia to the g-tensor. A (small) relativistic mass 
correction term AgfRMC [111 is included in both sets of data. 







LiiiGEir response 


This method 


A „( A M\ 


NeF 


Aoil 
II 


-336 


-328 


-414 




A 9± 


52633 


52778 


2935 


ArF 


A 3|| 


-349 


-343 


-4450 




Agx 


42439 


42519 


2914 


KrF 


A 9\\ 


-360 


-353 


-968 




Ag ± 


59920 


59674 


-1918 


XeF 


A 9\\ 


-358 


-354 


-3733 




A 9± 


163369 


158190 


-55099 


RnF 


A 9\\ 


-356 


-299 


-13670 




Ag ± 


603082 


488594 


-255079 



when benchmarking and/or developing improved DFT 
functionals for orbital magnetism. 

In the following we will show that the anisotropies in 
the orbital magnetizations are well described to allow 
us to calculate the electronic g-tensor of paramagnetic 
systems, in order to understand the microscopic structure 
of radicals or paramagnetic defects in solids. From the 
orbital magnetization we can obtain the deviation of the 
g-tensor, Ag M „ from the free electron value g e =2. 002319 
by the variation of M with a spin flip: 



A 9nv = -- 



M(e„) - M(-e„ 

s-(-s) 



ab 



M(e„) 
(14) 

where v, \i are Cartesian directions of the magnetic field 
and the total spin S, respectively. To get the full tensor 
Ag^ v , for every paramagnetic systems we carry out three 
calculations by aligning the spin quantization axis along 
the three Cartesian directions. 

To evaluate the approach, we compute the ^-tensors of 
selected diatomic radicals. An energy cutoff of 100 Ry is 
used in all molecular calculations. They are performed in 

TABLE III: Calculated principal values Ag in ppm for the 
diatomic molecules of the PbF-family. See Tab.HTlfor details. 







Linear Response 


This method 


Ag(AM) 


CF 


A 9\\ 


— DO 


-1999719 


-119746 




A 9± 


1920 


-553 


-240 


SiF 


A 9\\ 


— oo 


-1995202 


-100021 




Ag ± 


-480 


-2470 


-535 


GeF 


A 9\\ 


— oo 


-1998078 


-40609 




A 9± 


-15505 


-39101 


-388 


SnF 


A 9\\ 


— oo 


-1996561 


-72464 




Aff± 


-64997 


-142687 


-5339 


PbF 


A 9\\ 


— oo 


-1999244 


-90214 




A 9± 


-288383 


-556326 


-22476 
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a cubic repeated cell with a large volume of 8000 A 3 and 
the Brillouin zone is sampled only at the T point. For 
comparison, we also compute the ^-tensor via the lin- 
ear response method (LR) [Io| , which we recently imple- 
mented in the QuANTUM-Espresso package. For a wide 
range of molecular radicals including almost all of the 
examples discussed in Ref. flo| the new approach repro- 
duces the values obtained via LR within a few ppm (see 
also auxiliary Tab. I in (HJ). In Tab. ILTland lllll we report 
the calculated principal components of the computed ex- 
tensors for the RnF and PbF families. For the mem- 
bers of the RnF family qualitative deviations are only 
observed if heavy elements like Xe and Rn are involved, 
showing in LR large deviations Ag± of up to 10 5 ppm 
from g e for the corresponding fluorides. The treatment 
of SO-coupling beyond LR leads to considerably smaller 
values of Ag±, reduced by 3% (XeF) and 19% (RnF), re- 
spectively. Note that the reconstruction terms, Eqs. (TTTI 
IT3| . significantly contribute to the g-tensor. For the RnF 
family (see Tab. (TTJ) this is essential to obtain a value of 
Ag\\ ss 22 1 as also expected analytically (23l |. 

In contrast to the RnF family (5 electrons in the 



(Grant No. GE 1260/3-1) and by the CNRS. D. C. ac- 
knowledges partial support from ENI. Calculations were 
performed at the IDRIS, Paris (Grant No. 061202) and at 
CINECA, Bologna (Grant Supercalcolo 589046187069). 



p-shell, 



electronic configuration), the PbF family 



has only one electron within the p-shell. Without SO- 
coupling the unpaired electron occupies a degenerate e- 
level. Consequently, without SO, the HOMO-LUMO gap 
between the unpaired electron and the empty levels is 
zero, leading within LR to diverging values gu . This 
failure of LR is observed for all members of the PbF 
family, already for CF containing light elements exclu- 
sively. In contrast, our new method circumvents pertur- 
bation theory, and predicts a nearly vanishing g-value 
011 = 9e + Agii R3 along the bond direction of the di- 
atomic molecules as expected analytically (23| . 

In conclusion, we have shown how a recently devel- 
oped formula for the orbital magnetization can be ap- 
plied in an ab-initio pseudopotential scheme whereby the 
spin-orbit coupling enters explicitly the self-consistent cy- 
cle. In comparison with linear response methods, our ap- 
proach allows an improved calculation of the electronic 
g-tensor of paramagnetic systems containing heavy ele- 
ments or with large deviations of the g-tensor from the 
free electron value. The latter situation is encountered 
in many paramagnetic centers in solids, such as those 
exhibiting a Jahn- Teller distortion [24| and/or contain- 
ing transition metal impurities. In addition, our method 
provides improved orbital magnetizations with respect 
to the preexisting approaches that neglect the contribu- 
tions of the interstitial regions. This has been shown for 
the highly ordered ferromagnets where the orbital con- 
tribution is partially quenched by the crystal field. The 
presented approach is perfectly suited to describe also 
the ferromagnetism of nanostructures where the orbital 
quench is weaker and the orbital part of the magnetic 
moments becomes more dominant. 

U. G. acknowledges financial support by the DFG 
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In this auxiliary material we repo rt ( 1) the calculated EPR g-tensor for a set of molecular radicals including almost 
all of the examples discussed in Ref. [lOj; (2) the decomposition of the orbital magnetization of Fe, Co and Ni according 
to Eq. (9); (3) the orbital magnetization of Fe, Co and Ni, integrated within atomic spheres. 



radical 




LR 


This method 


g(M') 


g (AM) A ffRMC 




A <?|| 


-39.3 


-39.3 


24.7 





-64 




Ag± 


-41.7 


-41.7 


22.3 





-64 


CN 


A n 


-141 


-139 


32 


9 


-180 




A 9i 


-2600 


-2603 


-2192 


—231 


-180 


CO+ 


A n 


-136 


-134 


12 


33 


-179 






-3229 


-3231 


-3052 


—260 


-179 


BO 


A 9\\ 


-70 


-75 


-5 


22 


-92 




A 9± 


-2384 


-2384 


-2163 


— 129 


-92 


BS 


A.9|| 


-81 


-82 


-154 


177 


-105 






-9990 


-10001 


-9513 


—382 


-105 


4 in 


^9\\ 


1 /LQ 


-149 


339 


-294 


-192 




A 9x 


-1834 


-1842 


-1316 


-334 


-192 


NeF 


A 9\\ 


-336 


-328 


86 


6 


-420 




A.g± 


52633 


52778 


49843 


3355 


-420 


MgF 


A 9\\ 


-59 


-68 


57 


-37 


-88 




A 9x 


-2283 


-2316 


-2227 


-1 


-88 


ArF 


A<?|| 


-349 


-343 


102 


-10 


-435 




A 9x 


42439 


42519 


39605 


3349 


-435 


KrF 


A 9\\ 


-360 


-353 


615 


-520 


-448 




A 9± 


59920 


59674 


61593 


-1470 


-448 


XeF 


A 9\\ 


-358 


-354 


3380 


-3283 


-450 




A 9x 


163369 


158190 213285 


-54649 


-450 


HgF 


A 9\\ 


-288 


-263 


54601 


-54490 


-374 




A 9± 


-34268 


-33355 


52161 


-85115 


-374 


RnF 


A 9\\ 


-356 


-299 


13371 


-13196 


-474 




A 9± 


603082 


488594 743638 


-254605 


-474 



TABLE I: calculated At? in ppm for diatomic molecules, by linear response (LR) [lCf and with the current method. 
For sake of comparison, the SOO contribution is omitted from the GIPAW results. The "A contrib." column contains 
the sum of the contributions of AMbarc, AM para and AMdia to the g-tensor. The relativistic mass correction term 
A^rmc included in both sets of data is given explicitly. 
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Metal 




Mbare 


AAibarc 


A JVfpara 


AA/dia 


bcc-Fe 


LDA 


0.0616 


0.0005 


0.0016 


0.0003 




PBE 


0.0639 


0.0000 


0.0016 


0.0003 


fcc-Co 


LDA 


0.0706 


0.0019 


0.0014 


0.0002 




PBE 


0.0722 


0.0018 


0.0014 


0.0002 


hep-Co 


LDA 


0.0875 


0.0032 


0.0014 


0.0003 




PBE 


0.0908 


0.0032 


0.0014 


0.0003 


/cc-Ni 


LDA 


0.0519 


0.0019 


0.0007 


0.0000 




PBE 


0.0494 


0.0017 


0.0007 


0.0001 



TABLE II: Contributions to the orbital magnetization along the easy axis, in /zb per atom. See eq. (9) in the text. 
As in the case of molecules, the "A contrib." is not negligible and it is comparable to the difference between the full 
orbital magnetization and the orbital magnetization calculated inside atomic spheres (see Tab. Ill in this auxiliary 
material. 



Metal 


FLAPW LDA [5] FLAPW PBE [5] 


This work LDA This work PBE 


6cc-Fe 


0.048 


0.045 


0.0452 0.0433 


fcc-Co 


0.076 


0.073 


0.0641 0.0634 


hep-Co 






0.0835 0.0868 


/cc-Ni 


0.049 


0.050 


0.0499 0.0511 



TABLE III: Orbital magnetization contribution inside atomic spheres, in /ib per atom, along the easy axis. This is 
defined as M* rb = J n u* k (r) r x (— zV + k)u„k(r) dr where Q s is a sphere centered on one atom, of radius -Rmt- 
Rmt is given in units of the Bohr radius (ao). Rmt — 2.0 ao is a typical muffin-tin radius used in FLAPW calculations 
for Fe, Co and Ni. Our results agree very well with FLAPW results. By comparing to the orbital magnetization 
calculated according to the periodic formula (Tab. Ill of the paper), which takes into account not only the atomic 
spheres but also the interstitial region, it is evident that con tribution from the interstitial is not negligible. 



